from time import time
import numpy as np
import math
from pyXSteam.XSteam import XSteam
from scipy.integrate import solve_ivp



# %%  汽轮机模型
# 汽轮机机理模型
class Steam_turbine:
    def __init__(self):
        self.Xsteam = XSteam(XSteam.UNIT_SYSTEM_MKS)  # 97 公式计算
        # 模型参数(初始值)
        self.K_list = [12.3 / 513.6, 6.9 / 583, 4 / 550, 2 / 520, 0.9 / 500, 0.53 / 485, 0.21 / 450, 0.08 / 420]  # 流量系数列表
        self.ita_c = 0.90  # 汽轮机末级效率  (0-1)
        self.ita_t = 0.90  # 汽轮机的机械效率 (0-1)

    def update_parameters(self, ita_c, ita_t, k_list):
        """更新参数"""
        self.ita_c = ita_c
        self.ita_t = ita_t
        self.K_list = k_list

    """汽轮机-调节级"""

    def regulating_level(self, Tst, Pst, Dst, k_list):
        """ 根据主蒸汽流量 Dst 的不同范围，计算调节级后的蒸汽参数和内能
        :param Tst: 主蒸汽温度
        :param Pst: 主蒸汽压力 (MPa)
        :param Dst: 主蒸汽流量 (t/h)
        :param k_list: 流量系数列表（用于计算压力）
        :return: Dt(调节级后蒸汽流量); pt(调节级后压力); ht(调节级后焓值); t(调节级后温度); ne(内能)
        """
        k_reg = k_list[0]  # 调节级流量系数
        # 常数
        pt_vwo = 14  # VWO（全阀开启）工况下调节级后压力 (MPa)
        pstb_vwo = 16.67  # VWO工况下部分开启调节阀后的压力 (MPa)
        Dstb_vwo = 290.76  # 设计工况下部分开启调节阀的蒸汽流量 (t/h)
        # pt_d = 12.3      # 设计工况下调节级后压力
        # Dst_d = 513.6    # 设计工况下主蒸汽流量
        fai = 1  # 喷嘴和动叶速度系数（无单位）
        Kn = 1.3  # 蒸汽绝热指数（无单位）
        Dst_vwo = 581.5138889  # VWO工况下主蒸汽流量 (t/h)
        hst = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值

        """调节级:根据主蒸汽流量 Dst 分为三种情况"""
        # """(低流量, 仅部分开启调节阀)"""
        if Dst <= 290.7569444:
            Dt = Dst  # 调节级后流量等于主蒸汽流量
            pt = Dst * k_reg  # 调节级后压力与流量成正比
            pst2 = Dst * ((pstb_vwo - pt_vwo) / Dstb_vwo) + pt  # 蒸汽流经部分开调节阀后的压力(部分开启阀门后的压力，通过线性插值计算)
            hst2 = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值
            vst2 = self.Xsteam.v_ph(pst2 * 10, hst2)  # 部分开启阀门后的比体积 (m³/kg)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数，考虑绝热膨胀
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度(经验公式计算温度)
            ht = self.Xsteam.h_pt(pt * 10, tt2)  # 调节级后焓值
            t = self.Xsteam.t_ph(pt * 10, ht)  # 调节级后温度

        # """(低流量, 仅部分开启调节阀)"""
        elif 290.7569444 < Dst <= 436.1354167:
            # 部分开启阀门
            pt = Dst * k_reg  # 变工况下调节级后压力
            pst2 = (Dst - 0.5 * Dst_vwo) * ((pstb_vwo - pt_vwo) / Dstb_vwo) + pt  # 蒸汽流经部分开调节阀后的压力
            hst2 = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值
            vst2 = self.Xsteam.v_ph(pst2 * 10, hst2)  # 部分开启阀门后的比体积 (m³/kg)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数，考虑绝热膨胀
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度
            ht2 = self.Xsteam.h_pt(pt * 10, tt2)  # 通过部分开启调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值ht(部分开启阀门后焓值)
            # 全开阀门
            hst1 = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值
            pst1 = 0.95 * Pst  # 全开阀门后压力（假设 5% 压降）
            vst1 = self.Xsteam.v_ph(pst1 * 10, hst1)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数
            vt1 = (pst1 / pt) ** (1 / nm) * vst1  # 通过全开调节阀蒸汽流经喷嘴动叶后比体积
            tt1 = (1000 * pt * vt1 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过全开调节阀蒸汽流经喷嘴动叶后温度
            ht1 = self.Xsteam.h_pt(pt * 10, tt1)  # 通过全开调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值
            # 混合
            Dt = Dst  # 调节级后蒸汽流量
            ht = ((0.5 * Dst_vwo) * ht1 + (Dst - 0.5 * Dst_vwo) * ht2) / Dst  # 混合后调节级蒸汽焓值ht(按流量比例混合焓值)
            t = self.Xsteam.t_ph(pt * 10, ht)

            # 高流量(类似中流量，但混合比例调整为 75% 全开阀门流量)
        else:
            # 部分开启阀门
            pt = Dst * k_reg  # 变工况下调节级后压力
            pst2 = (Dst - 0.75 * Dst_vwo) * ((pstb_vwo - pt_vwo) / Dstb_vwo) + pt  # 蒸汽流经部分开调节阀后的压力
            hst2 = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值
            vst2 = self.Xsteam.v_ph(pst2 * 10, hst2)  # 流经部分开阀后的蒸汽的比体积
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数，考虑绝热膨胀
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度
            ht2 = self.Xsteam.h_pt(pt * 10, tt2)  # 通过部分开启调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值ht
            # 全开阀门
            hst1 = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值
            pst1 = 0.95 * Pst  # 蒸汽流经全开调节阀门后压力pst1
            vst1 = self.Xsteam.v_ph(pst1 * 10, hst1)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数
            vt1 = (pst1 / pt) ** (1 / nm) * vst1  # 通过全开调节阀蒸汽流经喷嘴动叶后比体积
            tt1 = (1000 * pt * vt1 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过全开调节阀蒸汽流经喷嘴动叶后温度
            ht1 = self.Xsteam.h_pt(pt * 10, tt1)  # 通过全开调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值
            # 混合
            Dt = Dst  # 调节级后蒸汽流量
            ht = ((0.75 * Dst_vwo) * ht1 + (Dst - 0.75 * Dst_vwo) * ht2) / Dst  # 混合后调节级蒸汽焓值ht
            t = self.Xsteam.t_ph(pt * 10, ht)
        # 内能 (kW)，表示调节级能量转换。
        ne = Dst * (hst - ht)
        return Dt, pt, ht, t, ne

    """汽轮机-压力级"""
    "压力级第一级"

    def pressure_level_1(self, Dt, pt, ht, k_list):
        """
        :param Dt: 入口蒸汽流量 (t/h)
        :param pt: 入口压力 (MPa)
        :param ht: 入口焓值 (kJ/kg)
        :param k_list: 流量系数列表
        :return: D1(抽汽流量); p1(抽汽压力); h1(抽汽焓值); ne(内能)
        """
        k_pre_1 = k_list[1]
        p1d = 6.90877688  # 设计工况抽汽压力
        xl = 0.9  # 级间效率
        Dj1d = 583.3333333  # 设计工况下第一级蒸汽流量
        p1 = Dt * k_pre_1  # 第一级抽汽压力，即第一级出口蒸汽压力(抽汽压力与流量成正比)
        s1 = self.Xsteam.s_ph(pt * 10, ht)  # 级组入口蒸汽熵
        h11 = self.Xsteam.h_ps(p1 * 10, s1)  # 第一级抽汽理想焓值(等熵膨胀)
        h1 = ht - xl * (ht - h11)  # 第一级实际抽汽焓值(实际焓值，考虑效率)
        Pa = 0.101  # 大气压 (MPa)
        EGJCQ = 3
        ps1 = (p1 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps1 > p1:
            D1 = 0
        else:
            D1 = math.sqrt(p1 - ps1) * (-0.4972 * p1 ** 3 + 7.227 * p1 ** 2 - 19.44 * p1 + 48.69)  # 第一级抽汽流量
        ne = Dt * (ht - h1)  # 内能
        return D1, p1, h1, ne

    "压力级第二级"

    def pressure_level_2(self, D1, p1, h1, Dt, k_list):
        k_pre_2 = k_list[2]
        p2d = 4.065  # 设计工况抽汽压力
        xl = 0.9  # 级间效率
        Dj2d = 550.548348  # 设计工况下第二级蒸汽流量
        Dj2 = Dt - D1  # 第二级蒸汽流量
        p2 = Dj2 * k_pre_2  # 第二级抽汽压力，即第二级出口蒸汽压力
        s2 = self.Xsteam.s_ph(p1 * 10, h1)  # 级组二入口蒸汽熵
        h22 = self.Xsteam.h_ps(p2 * 10, s2)  # 第二级抽汽理想焓值
        h2 = h1 - xl * (h1 - h22)  # 第二级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps2 = (p2 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps2 > p2:
            D2 = 0
        else:
            D2 = math.sqrt(p2 - ps2) * (-2.75 * p2 ** 3 + 23.6 * p2 ** 2 - 38.75 * p2 + 62.27)
        ne = Dj2 * (h1 - h2)
        return D2, p2, h2, Dj2, ne

    "压力级第三级"

    def pressure_level_3(self, Pre, Tre, Dre, k_list):
        """
        :param Pre: 再热蒸汽压力
        :param Tre: 再热蒸汽温度
        :param Dre: 再热蒸汽流量
        :param k_list: 流量系数
        :return:
        """
        k_pre_3 = k_list[3]
        p3d = 2.090  # 设计工况下第三级抽汽压力
        xl = 0.9  # 级间效率
        Dj3d = 520.07  # 设计工况下第三级蒸汽流量
        p3 = Dre * k_pre_3  # 第三级抽汽压力，即第三级出口蒸汽压力
        hre = self.Xsteam.h_pt(Pre * 10, Tre)  # 级组三入口蒸汽焓值
        s3 = self.Xsteam.s_ph(Pre * 10, hre)  # 级组三入口蒸汽熵值
        h33 = self.Xsteam.h_ps(p3 * 10, s3)  # 第三级抽汽理想焓值
        h3 = hre - xl * (hre - h33)  # 第三级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps3 = (p3 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps3 > p3:
            D3 = 0
        else:
            D3 = math.sqrt(p3 - ps3) * (-6.841 * p3 ** 3 + 31.39 * p3 ** 2 - 17.2 * p3 + 22.72)  # 第三级抽汽流量
        Dj3 = Dre - D3  # 第三级蒸汽流量
        ne = Dre * (hre - h3)
        return D3, p3, h3, Dj3, ne, hre

    "压力级第四级"

    def pressure_level_4(self, h3, p3, Dj3, k_list):
        k_pre_4 = k_list[4]
        p4d = 0.916990  # 设计工况下第四级抽汽压力
        xl = 0.9  # 级间效率
        Dj4d = 500.07  # 设计工况下第四级蒸汽流量
        p4 = Dj3 * k_pre_4  # 第四级抽汽压力，即第四级出口蒸汽压力
        s4 = self.Xsteam.s_ph(p3 * 10, h3)  # 级组四入口蒸汽熵
        h44 = self.Xsteam.h_ps(p4 * 10, s4)  # 第四级抽汽理想焓值
        h4 = h3 - xl * (h3 - h44)  # 第四级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps4 = (p4 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps4 > p4:
            D4 = 0
        else:
            D4 = math.sqrt(p4 - ps4) * (-64.69 * p4 ** 3 + 144.7 * p4 ** 2 - 47.91 * p4 + 35.01)  # 第四级抽汽流量
        ne = Dj3 * (h3 - h4)
        Dj4 = Dj3 - D4  # 第四级蒸汽流量

        return D4, p4, h4, Dj4, ne

    "压力级第五级"

    def pressure_level_5(self, h4, p4, Dj4, k_list):
        k_pre_5 = k_list[5]
        p5d = 0.53591  # 设计工况下第五级抽汽压力
        xl = 0.9  # 级间效率
        Dj5d = 485.54  # 设计工况下第五级蒸汽流量
        p5 = Dj4 * k_pre_5  # 第五级抽汽压力，即第五级出口蒸汽压力
        s5 = self.Xsteam.s_ph(p4 * 10, h4)  # 级组五入口蒸汽熵
        h55 = self.Xsteam.h_ps(p5 * 10, s5)  # 第五级抽汽理想焓值
        h5 = h4 - xl * (h4 - h55)  # 第五级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps5 = (p5 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps5 > p5:
            D5 = 0
        else:
            D5 = math.sqrt(p5 - ps5) * (-365.1 * p5 ** 3 + 529.1 * p5 ** 2 - 68.05 * p5 + 56.83)  # 第五级抽汽流量
        ne = Dj4 * (h4 - h5)
        Dj5 = Dj4 - D5  # 第五级蒸汽流量

        return D5, p5, h5, Dj5, ne

    "压力级第六级"

    def pressure_level_6(self, h5, p5, Dj5, k_list):
        k_pre_6 = k_list[6]
        p6d = 0.21902712  # 设计工况下第六级抽汽压力
        xl = 0.9  # 级间效率
        Dj6d = 450.07  # 设计工况下第六级蒸汽流量
        p6 = Dj5 * k_pre_6  # 第六级抽汽压力，即第六级出口蒸汽压力
        s6 = self.Xsteam.s_ph(p5 * 10, h5)  # 级组六入口蒸汽熵
        h66 = self.Xsteam.h_ps(p6 * 10, s6)  # 第六级抽汽理想焓值
        h6 = h5 - xl * (h5 - h66) + 110  # 第六级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps6 = (p6 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps6 > p6:
            D6 = 0
        else:
            D6 = math.sqrt(p6 - ps6) * (9233 * p6 ** 3 - 4349 * p6 ** 2 + 1313 * p6 + 2.72)  # 第六级抽汽流量
        ne = Dj5 * (h5 - h6)
        Dj6 = Dj5 - D6  # 第六级蒸汽流量
        return D6, p6, h6, Dj6, ne

    "压力级第七级"

    def pressure_level_7(self, h6, p6, Dj6, k_list):
        k_pre_7 = k_list[7]
        p7d = 0.081204  # 设计工况下第七级抽汽压力
        xl = 0.9  # 级间效率
        Dj7d = 420.95  # 设计工况下第七级蒸汽流量
        p7 = Dj6 * k_pre_7  # 第七级抽汽压力，即第七级出口蒸汽压力
        s7 = self.Xsteam.s_ph(p6 * 10, h6)  # 级组七入口蒸汽熵
        h77 = self.Xsteam.h_ps(p7 * 10, s7)  # 第七级抽汽理想焓值
        h7 = h6 - xl * (h6 - h77)  # 第七级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps7 = (p7 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps7 > p7:
            D7 = 0
        else:
            D7 = math.sqrt(p7 - ps7) * (5608 * p7 ** 3 - 604 * p7 ** 2 + 2841 * p7 - 54.62)  # 第七级抽汽流量
        ne = Dj6 * (h6 - h7)
        Dj7 = Dj6 - D7  # 第七级蒸汽流量
        return D7, p7, h7, Dj7, ne

    "末级"

    def final_level(self, h7, p7, Dj7, Pc):
        """
        :param h7: 压力级出口焓值（第七级焓值）
        :param p7: 压力级出口压力（第七级压力）
        :param Dj7: 压力级流量（第七级流量）
        :param Pc: 背压
        :return: hc（排汽焓值）; Dc (排汽流量); ne(内能); Tp(排汽温度); Tl_p(饱和温度)
        """
        # 汽轮机末级效率
        xl = self.ita_c
        Dc = Dj7  # 末级蒸汽流量
        sm = self.Xsteam.s_ph(p7 * 10, h7)  # 末级入口蒸汽熵
        hm = self.Xsteam.h_ps(Pc * 10, sm)  # 末级抽汽理想焓值
        hc = h7 - xl * (h7 - hm)  # 末级排汽焓值
        ne = Dc * (h7 - hc)  # 内能
        Tc = self.Xsteam.t_ph(Pc * 10, hc)  # 排汽温度

        return hc, Dc, Tc, ne

    "汽轮机仿真：综合模拟汽轮机各阶段，计算负荷"

    def Turbine_simulation(self, Tst, Pst, Dst, Pre, Tre, Pc):
        """
        :param Tst: 主蒸汽温度
        :param Pst: 主蒸汽压力
        :param Dst: 主蒸汽流量
        :param Pre: 再热蒸汽压力
        :param Tre: 再热蒸汽温度
        :param Pc: 背压
        :return: Dt(调节级后流量); Pt(调节级压力); Dc(末级排汽流量); Ne(总负荷)
        """
        k_list = self.K_list  # 流量系数列表
        # 调节级
        Dt, Pt, ht, t, ne_regulate = self.regulating_level(Tst, Pst, Dst, k_list)
        # 压力级第一级
        D1, p1, h1, ne1 = self.pressure_level_1(Dt, Pt, ht, k_list)
        # 压力级第二级
        D2, p2, h2, Dj2, ne2 = self.pressure_level_2(D1, p1, h1, Dt, k_list)
        # 压力级第三级 (再热器)
        D3, p3, h3, Dj3, ne3, hre = self.pressure_level_3(Pre, Tre, Dj2, k_list)
        # 压力级第四级
        D4, p4, h4, Dj4, ne4 = self.pressure_level_4(h3, p3, Dj3, k_list)
        # 压力级第五级
        D5, p5, h5, Dj5, ne5 = self.pressure_level_5(h4, p4, Dj4, k_list)
        # 压力级第六级
        D6, p6, h6, Dj6, ne6 = self.pressure_level_6(h5, p5, Dj5, k_list)
        # 压力级第七级
        D7, p7, h7, Dj7, ne7 = self.pressure_level_7(h6, p6, Dj6, k_list)
        # 末级
        hc, Dc, Tc, ne_final = self.final_level(h7, p7, Dj7, Pc)
        # 总负荷 (总内能乘以机械效率并转换为 MW)
        Ne = (ne_regulate + ne1 + ne2 + ne3 + ne4 + ne5 + ne6 + ne7 + ne_final) * self.ita_t
        Ne = Ne * 0.001  # 转换为 MW

        return Dt, Pt, Dc, hc, Tc, Ne


# %% 直接空冷凝汽器的机理模型
class Air_Condenser:
    def __init__(self):
        self.Xsteam = XSteam(XSteam.UNIT_SYSTEM_MKS)
        # 凝汽器的辨识参数
        self.kc = 0.04  # 凝汽器换热系数
        self.cf = 0.9  # 风机集群因子
        self.kec = 1.5  # 凝结水流量修正系数

        # 凝汽器的模型参数
        self.Vn = 2.5  # 额定风速 (m/s)
        self.Ac = 1440216  # 换热表面积 (m²)
        self.Af = 10731  # 风机覆盖面积 (m²)
        self.rho_air = 1.043  # 空气密度 (kg/m³)
        self.Cp_air = 1.005  # 空气比热容 (kJ/(kg·℃))
        self.Qvd = 550  # 风机设计体积流量 (m³/s)
        self.n_design = 110  # 设计转速
        self.Rs = 461.5  # 水蒸气气体常数 (J/(kg·K))
        self.Vc = 9600 * 56  # 凝汽器体积 (m³)

    def update_parameters(self, kc, kec):
        """更新参数"""
        self.kc = kc
        self.kec = kec

    def Air_Condenser_Sim(self, t, x, u):
        """定义空冷凝汽器的微分方程，用于 solve_ivp 求解"""
        Pc = x[0]  # 背压 (MPa)
        Dc, hc, Ta, fan_speed, fan_number = u  # 控制输入：排汽流量、排汽焓值、环境温度、风机转速、风机运行台数

        k = self.kc
        cluster_factor = self.cf

        # 凝汽器的模型参数
        Vn = self.Vn
        Ac = self.Ac
        Af = self.Af
        rho_air = self.rho_air
        Cp_air = self.Cp_air
        Qvd = self.Qvd
        n_design = self.n_design
        Rs = self.Rs
        Vc = self.Vc
        kec = self.kec

        pc = Pc * 10  # MPa to bar
        hc_sa = self.Xsteam.hV_p(pc)  # 凝汽器压力对应的饱和蒸汽焓(理想排汽焓值)
        Tc_sa = self.Xsteam.tsat_p(pc)  # 凝汽器压力对应的饱和温度
        hlc = self.Xsteam.hL_p(pc)  # 凝汽器压力对应的饱和水焓

        # 计算风机转速
        Qv = Qvd * (fan_speed / n_design)  # 单台风机体积流量
        Vf_current = (Qv * fan_number * cluster_factor) / Af  # 当前风速
        Vf = Vf_current

        # kc = k * (np.power(complex(Vf / Vn), 0.6)).real  # 修正后的传热系数
        NTU = (k * Ac) / (rho_air * Vf * Af * Cp_air)  # 传热单元数NTU
        Qc = (rho_air * Vf * Af * Cp_air) * (1 - math.exp(-NTU)) * (Tc_sa - Ta)  # 此汽轮机排汽和环境空气之间的换热量 Qc
        Dec = Qc / (hc - hlc) * kec  # 凝汽器内凝结的排汽质量流量
        dPc_dt = (Rs * (Tc_sa + 273.15) * (Dc - Dec)) / Vc  # Pc 的更新量（单位为：Pa）
        dPc_dt = dPc_dt * 0.000001  # Pa → MPa

        return np.array([dPc_dt])

    def Calculate_current_state(self, Pc, Dc, hc, Ta, fan_speed, fan_number):
        """计算当前时刻的凝结水流量和风速"""
        pc = Pc * 10  # MPa to bar
        hc_sa = self.Xsteam.hV_p(pc)  # 凝汽器压力对应的饱和蒸汽焓(理想排汽焓值)
        Tc_sa = self.Xsteam.tsat_p(pc)  # 凝汽器压力对应的饱和温度
        hlc = self.Xsteam.hL_p(pc)  # 凝汽器压力对应的饱和水焓

        k = self.kc
        cluster_factor = self.cf
        Vn = self.Vn
        Ac = self.Ac
        Af = self.Af
        rho_air = self.rho_air
        Cp_air = self.Cp_air
        Qvd = self.Qvd
        n_design = self.n_design
        kec = self.kec

        # 计算风机转速
        Qv = Qvd * (fan_speed / n_design)  # 单台风机体积流量
        Vf_current = (Qv * fan_number * cluster_factor) / Af  # 当前风速
        Vf = Vf_current

        kc = k * (np.power(complex(Vf / Vn), 0.6)).real  # 修正后的传热系数
        NTU = (kc * Ac) / (rho_air * Vf * Af * Cp_air)  # 传热单元数NTU
        Qc = (rho_air * Vf * Af * Cp_air) * (1 - math.exp(-NTU)) * (Tc_sa - Ta)  # 此汽轮机排汽和环境空气之间的换热量 Qc
        Dec = Qc / (hc - hlc) * kec  # 凝汽器内凝结的排汽质量流量

        return Vf, Dec, Pc, kc

    def Condenser_Simulation(self, Pc, Dc, hc, Ta, fan_speed, fan_number, t_end=60.0, dt=0.1):
        """ 执行一步积分仿真，返回当前时刻参数以及积分后下一时刻的背压
            t_end: 仿真时间段长度 (默认60秒，即采样时间)
            dt: 积分步长 (默认0.1秒)
        """
        # 计算当前时刻的状态(风速；凝结水流量；背压；kc参数)
        Vf, Dec, Pc, kc = self.Calculate_current_state(Pc, Dc, hc, Ta, fan_speed, fan_number)

        tspan = np.arange(0, t_end + dt, dt)
        u_ctrl = np.array([Dc, hc, Ta, fan_speed, fan_number])
        state = np.array([Pc])

        # 背压模型进行积分
        sol = solve_ivp(
            fun=lambda t, x: self.Air_Condenser_Sim(t, x, u_ctrl),
            t_span=(tspan[0], tspan[-1]),
            y0=state,
            t_eval=tspan)
        Pc_next = sol.y[:, -1]  # 取最后一个作为下一时刻的状态
        Pc_next = Pc_next[0]

        return Pc_next, Dec, Pc


# %% 风机功耗模型
class Air_Power:
    def __init__(self):
        pass

    @staticmethod
    def Air_Power_Sim(fan_speed, fan_number):
        # 计算风机功耗
        n0 = 110  # 额定风机转速
        P0 = 115  # 额定风机功率(115KW)
        Power_fan = ((fan_speed / n0) ** 3) * P0 * fan_number * 0.001  # 单位:MW
        return Power_fan


#%% 耦合建模
class Model_Simulation:
    def __init__(self, config):
        """ 初始化模拟类，汽轮机仿真 """
        self.Xsteam = XSteam(XSteam.UNIT_SYSTEM_MKS)
        self.Turbine = Steam_turbine()  # 蒸汽轮机模型
        self.Condenser = Air_Condenser()  # 空冷凝汽器模型
        self.Condenser_Power = Air_Power()  # 风机功耗模型

        # 在线梯度下降更新效率参数，滑动窗口最小二乘法更新流量系数
        self.ita_c = self.Turbine.ita_c
        self.ita_t = self.Turbine.ita_t
        self.k_list = self.Turbine.K_list
        self.learning_rate = 2.5e-6  # 学习率
        self.P_buf = []
        self.D_buf = []
        self.kc = self.Condenser.kc
        self.kec = self.Condenser.kec

        self.Pc_current = None  # 初始化当前时刻背压
        self.mech_pre = None
        self.Pc_pre = None

        # 步骤计数器和重置间隔
        self.step_count = 0
        self.reset_interval = 300  # 每隔300步用真实值替换 Pc_current

    def Sim_Step(self, data, Pc):
        """ 执行单步模拟仿真 """
        Tst, Pst, Dst, Pre, Tre, Ta, fan_speed, fan_number = data
        Dt, Pt, Dc, hc, Tc, Ne = self.Turbine.Turbine_simulation(Tst, Pst, Dst, Pre, Tre, Pc)
        Pc_next, Dec, Pc = self.Condenser.Condenser_Simulation(Pc, Dc, hc, Ta, fan_speed, fan_number)

        return np.array([Dt, Pt, Dc, hc, Tc, Ne, Pc, Dec, Pc_next])

    @staticmethod
    # 根据抽汽压力计算实际抽汽流量(工程计算方法)
    def Di_calc(p_real):
        """ 计算每个压力级的抽汽流量 """
        Pa = 0.101  # 大气压
        EGJCQ = 3  # 抽汽损失百分比, 3 %
        pressures = p_real
        Di_real = []
        for p in pressures:
            ps_real = (p + Pa) * (1 - EGJCQ / 100) - Pa
            D_real = math.sqrt(p - ps_real) * (-0.4972 * p ** 3 + 7.227 * p ** 2 - 19.44 * p + 48.69)
            Di_real.append(D_real)
        return tuple(Di_real)

    @staticmethod
    # 汽轮机流量系数计算(滑动窗口最小二乘法)
    def turbine_coff(D, P):
        k = P @ np.linalg.pinv(D)
        return k

    def Update_Turbine_Parameters(self, mech_in, p_real, Ne_k, Pc_k):
        """  使用梯度下降更新模型参数 ita_c 和 ita_t 以及流量系数列表
        :param mech_in: 当前输入
        :param p_real: 当前实际抽气压力(调节级+7个压力级)
        :param Ne_k: 当前时刻负荷
        :param Pc_k: 当前时刻背压
        :return:
        """
        # 使用当前参数运行模型
        self.Turbine.update_parameters(self.ita_c, self.ita_t, self.k_list)
        self.Condenser.update_parameters(self.kc, self.kec)

        states = self.Sim_Step(mech_in, Pc_k)
        Ne_sim = states[5]  # 仿真的机组负荷
        error_Ne = Ne_sim - Ne_k  # Ne误差

        # 数值梯度估计(有限差分法)
        delta_ita_c = self.ita_c / 1000
        delta_ita_t = self.ita_t / 1000

        # 计算 ita_c 的梯度
        ita_c_plus = self.ita_c + delta_ita_c
        self.Turbine.update_parameters(ita_c_plus, self.ita_t, self.k_list)
        states_plus_ita_c = self.Sim_Step(mech_in, Pc_k)
        Ne_plus = states_plus_ita_c[5]
        ita_c_minus = self.ita_c - delta_ita_c
        self.Turbine.update_parameters(ita_c_minus, self.ita_t, self.k_list)
        states_minus_ita_c = self.Sim_Step(mech_in, Pc_k)
        Ne_minus = states_minus_ita_c[5]
        gradient_ita_c = (Ne_plus - Ne_minus) / (2 * delta_ita_c)

        # 计算 ita_t 的梯度
        ita_t_plus = self.ita_t + delta_ita_t
        self.Turbine.update_parameters(self.ita_c, ita_t_plus, self.k_list)
        states_plus_ita_t = self.Sim_Step(mech_in, Pc_k)
        Ne_plus = states_plus_ita_t[5]
        ita_t_minus = self.ita_t - delta_ita_t
        self.Turbine.update_parameters(self.ita_c, ita_t_minus, self.k_list)
        states_minus_ita_t = self.Sim_Step(mech_in, Pc_k)
        Ne_minus = states_minus_ita_t[5]
        gradient_ita_t = (Ne_plus - Ne_minus) / (2 * delta_ita_t)

        # 更新参数
        self.ita_c -= self.learning_rate * error_Ne * gradient_ita_c
        self.ita_c = max(0.75, min(1.0, self.ita_c))  # 限制范围

        self.ita_t -= self.learning_rate * error_Ne * gradient_ita_t
        self.ita_t = max(0.75, min(1.0, self.ita_t))  # 限制范围

        # [Dt, Pt, Dc, hc, Tc, Ne, Pc, Dec, Pc_next]
        Dt_sim, Pt_sim, Dc_sim, hc_sim, Tc_sim, Ne_sim, Pc_sim, Dec_sim, Pc_next_sim = states
        Dj1_cur = Dt_sim
        # 计算压力级的实际抽汽流量
        D1_real, D2_real, D3_real, D4_real, D5_real, D6_real, D7_real = self.Di_calc(p_real[1:])
        # 实际的抽汽压力（调节级+7个压力级）
        P_real = p_real
        # 计算级间流量
        Dj2_real = Dj1_cur - D1_real
        Dj3_real = Dj2_real - D2_real
        Dj4_real = Dj3_real - D3_real
        Dj5_real = Dj4_real - D4_real
        Dj6_real = Dj5_real - D5_real
        Dj7_real = Dj6_real - D6_real
        D_real = np.asarray([[Dt_sim, Dj1_cur, Dj2_real, Dj3_real, Dj4_real, Dj5_real, Dj6_real, Dj7_real]])

        # 添加到缓冲区(滑动窗口设置为30min)
        if len(self.P_buf) < 30:
            self.P_buf.append(P_real)
            self.D_buf.append(D_real)
        else:
            # 移除最旧的数据
            self.P_buf.pop(0)
            self.D_buf.pop(0)
            self.P_buf.append(P_real)
            self.D_buf.append(D_real)
            # 将缓冲区转换为 NumPy 数组
            P_buf_array = np.vstack(self.P_buf)
            D_buf_array = np.vstack(self.D_buf)

            # 流量系数校正
            k0 = self.turbine_coff(D_buf_array[:, 0].reshape(1, -1), P_buf_array[:, 0].reshape(1, -1))
            k1 = self.turbine_coff(D_buf_array[:, 1].reshape(1, -1), P_buf_array[:, 1].reshape(1, -1))
            k2 = self.turbine_coff(D_buf_array[:, 2].reshape(1, -1), P_buf_array[:, 2].reshape(1, -1))
            k3 = self.turbine_coff(D_buf_array[:, 3].reshape(1, -1), P_buf_array[:, 3].reshape(1, -1))
            k4 = self.turbine_coff(D_buf_array[:, 4].reshape(1, -1), P_buf_array[:, 4].reshape(1, -1))
            k5 = self.turbine_coff(D_buf_array[:, 5].reshape(1, -1), P_buf_array[:, 5].reshape(1, -1))
            k6 = self.turbine_coff(D_buf_array[:, 6].reshape(1, -1), P_buf_array[:, 6].reshape(1, -1))
            k7 = self.turbine_coff(D_buf_array[:, 7].reshape(1, -1), P_buf_array[:, 7].reshape(1, -1))
            self.k_list = np.array([k0, k1, k2, k3, k4, k5, k6, k7]).reshape(-1)

    def Update_Condenser_Parameters(self, mech_pre, Pc_pre, Pc_cur):
        """
        :param mech_pre: 上一时刻的输入
        :param Pc_pre: 上一时刻的背压
        :param Pc_cur: 当前时刻的背压 (kPa)
        :return: 参数更新
        """
        # 使用当前参数运行模型
        self.Turbine.update_parameters(self.ita_c, self.ita_t, self.k_list)
        self.Condenser.update_parameters(self.kc, self.kec)

        mech_out = self.Sim_Step(mech_pre, Pc_pre)
        Pc_cur_pre = mech_out[-1] * 1000  # kPa
        error_Pc = Pc_cur_pre - Pc_cur  # Pc误差

        # 计算 kc 的梯度
        delta_kc = self.kc / 1000
        kc_plus = self.kc + delta_kc
        self.Condenser.update_parameters(kc_plus, self.kec)
        states_plus_kc = self.Sim_Step(mech_pre, Pc_pre)
        Pc_plus = states_plus_kc[-1]
        kc_minus = self.kc - delta_kc
        self.Condenser.update_parameters(kc_minus, self.kec)
        states_minus_kc = self.Sim_Step(mech_pre, Pc_pre)
        Pc_minus = states_minus_kc[-1]
        gradient_kc = (Pc_plus - Pc_minus) / (2 * delta_kc)

        # 计算kec的梯度
        delta_kec = self.kec / 1000
        kec_plus = self.kec + delta_kec
        self.Condenser.update_parameters(self.kc, kec_plus)
        states_plus_kec = self.Sim_Step(mech_pre, Pc_pre)
        Pc_plus = states_plus_kec[-1]
        kec_minus = self.kec - delta_kec
        self.Condenser.update_parameters(self.kc, kec_minus)
        states_minus_kec = self.Sim_Step(mech_pre, Pc_pre)
        Pc_minus = states_minus_kec[-1]
        gradient_kec = (Pc_plus - Pc_minus) / (2 * delta_kec)

        # 更新参数
        self.kc -= self.learning_rate * error_Pc * gradient_kc
        self.kc = max(0.0, min(0.04, self.kc))  # 限制范围
        self.kec = self.kec - 4.0 * error_Pc * gradient_kec
        self.kec = max(0.5, min(5.0, self.kec))  # 限制范围

    @staticmethod
    def calculate_fan_info(input_dict):
        """计算风机运行台数和平均风机转速"""
        # 生成所有转速键
        speed_keys = [f'HD_HA_30MAG{i}0AN00{j}CS' for i in range(1, 9) for j in range(1, 8)]
        # 获取所有转速值，默认0如果不存在
        speeds = [input_dict.get(key, 0) for key in speed_keys]
        # 过滤运行的风机转速 (>1 rpm)
        running_speeds = [s for s in speeds if s > 1]
        fan_number = len(running_speeds)
        fan_speed = sum(running_speeds) / fan_number if fan_number > 0 else 0
        return fan_number, fan_speed

    def process_step(self, input_dict):
        """
        处理当前时刻的数据，返回仿真输出。
        :param input_dict: 当前时刻的输入字典。
        :return: output_dict: 当前时刻的仿真输出字典。
        """
        # 提取输入
        Tst = input_dict.get('主汽温度', 1.0)
        Pst = input_dict.get('调门开度', 1.0)
        Dst = input_dict.get('主汽流量', 1.0)
        Pre = input_dict.get('再热蒸汽压力', 1.0)
        Tre = input_dict.get('再热蒸汽温度', 1.0)
        Ne_k = input_dict.get('机组负荷', 1.0)
        Pc_real = input_dict.get('背压', 1.0)           # kPa
        Ta = input_dict.get('空冷大气温度', 1.0)
        Pt = input_dict.get('调节级压力', 1.0)
        P_one = input_dict.get('一段抽汽压力', 1.0)
        P_two = input_dict.get('二段抽汽压力', 1.0)
        P_three = input_dict.get('三段抽汽压力', 1.0)
        P_four = input_dict.get('四段抽汽压力', 1.0)
        P_five = input_dict.get('五段抽汽压力', 1.0)
        P_six = input_dict.get('六段抽汽压力', 1.0)
        P_seven = input_dict.get('七段抽汽压力', 1.0)
        
        # 计算风机运行台数和平均转速
        fan_number, fan_speed = self.calculate_fan_info(input_dict)
        mech_in = [Tst, Pst, Dst, Pre, Tre, Ta, fan_speed, fan_number]
        p_real = [Pt, P_one, P_two, P_three, P_four, P_five, P_six, P_seven]

        # 初始化 Pc_current 使用真实 Pc
        if self.Pc_current is None:
            self.Pc_current = Pc_real / 1000  # MPa
            self.Pc_pre = self.Pc_current
        # 每隔 reset_interval 步，用真实 Pc 替换 Pc_current
        self.step_count += 1             # 步骤计数
        if self.step_count % self.reset_interval == 0:
            self.Pc_current = Pc_real / 1000    # MPa

        # 使用当前参数运行模型
        self.Turbine.update_parameters(self.ita_c, self.ita_t, self.k_list)
        self.Condenser.update_parameters(self.kc, self.kec)

        # 单步仿真
        states = self.Sim_Step(mech_in, self.Pc_current)
        Fan_power_current = self.Condenser_Power.Air_Power_Sim(fan_speed, fan_number)

        Pc_next = states[-1]

        # 更新模型参数
        if self.mech_pre is not None:
            self.Update_Turbine_Parameters(mech_in, p_real, Ne_k, self.Pc_current)
            self.Update_Condenser_Parameters(self.mech_pre, self.Pc_pre, Pc_real)       # Pc_real in kPa

        # 更新 pre
        self.mech_pre = mech_in
        self.Pc_pre = self.Pc_current
        self.Pc_current = Pc_next

        # 输出字典
        output_dict = {
            '机组排汽流量': states[2],                      # 机组排汽流量仿真
            '机组排汽焓值': states[3],                  # 机组排汽焓值仿真
            '机组汽轮机效率': self.ita_t,                    # 机组汽轮机效率
            '机组冷端总电耗': Fan_power_current,  # 机组冷端总电耗
            '机组模型背压仿真误差': Pc_real - states[6] * 1000,        # 机组模型背压仿真误差
            '机组模型负荷仿真误差': Ne_k - states[5],                          # 机组模型负荷仿真误差
            '机组平均风机转速': fan_speed                        # 机组平均风机转速
        }

        return output_dict



    def mockData(self, input_dict):
        # 模拟一点值，适当加点随机数
        import random
        import time
        time.sleep(5)
        print("接收到输入数据")
        print(input_dict)
        output_dict = {
            '机组排汽流量': 500.0 + random.uniform(-10, 10),            # 模拟值
            '机组排汽焓值': 3200.0 + random.uniform(-30, 30),           # 模拟值
            '机组汽轮机效率': 0.88 + random.uniform(-0.01, 0.01),        # 模拟值
            '机组冷端总电耗': 1200.0 + random.uniform(-100, 100),        # 模拟值
            '机组模型背压仿真误差': random.uniform(-1, 1),              # 模拟值
            '机组模型负荷仿真误差': random.uniform(-5, 5),               # 模拟值
            '机组平均风机转速': 105.0 + random.uniform(-5, 5)            # 模拟值
        }
        print("计算完成")
        return output_dict


def _build_mock_input(step, model):
    """构造一份模拟数据输入。所有单位与模型期望一致。"""
    Tst = 540.0
    Pst = 16.5
    Dst = 500.0
    Pre = 3.6
    Tre = 540.0
    Ne_k = 300.0
    Pc_real = 10.0
    Ta = 25.0 + 2.0 * math.sin(0.1 * step)

    k_list = model.k_list
    Pt = Dst * float(k_list[0])
    P_one = 6.0
    P_two = 4.0
    P_three = 2.2
    P_four = 0.95
    P_five = 0.55
    P_six = 0.22
    P_seven = 0.08

    input_dict = {
        '主汽温度': Tst,
        '调门开度': Pst,
        '主汽流量': Dst,
        '再热蒸汽压力': Pre,
        '再热蒸汽温度': Tre,
        '机组负荷': Ne_k,
        '背压': Pc_real,
        '空冷大气温度': Ta,
        '调节级压力': Pt,
        '一段抽汽压力': P_one,
        '二段抽汽压力': P_two,
        '三段抽汽压力': P_three,
        '四段抽汽压力': P_four,
        '五段抽汽压力': P_five,
        '六段抽汽压力': P_six,
        '七段抽汽压力': P_seven,
    }

    base_speeds = [100.0, 105.0, 98.0, 110.0, 95.0, 102.0]
    assign_list = []
    for i in range(1, 3):
        for j in range(1, 4):
            assign_list.append((i, j))
    for idx, (i, j) in enumerate(assign_list[:len(base_speeds)]):
        key = f'HD_HA_30MAG{i}0AN00{j}CS'
        input_dict[key] = base_speeds[idx] + 1.5 * math.sin(0.2 * step + idx)

    return input_dict


def main():
    """使用模拟数据运行若干步，打印关键输出。"""
    config = {
        'learning_rate': 0.01,
        'reset_interval': 300,
        'ita_c': 0.8,
        'ita_t': 0.8,
        'kc': 0.02,
        'kec': 0.5,
        'k_list': [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08],
    }
    model = Model_Simulation(config)
    steps = 5
    for step in range(1, steps + 1):
        inp = _build_mock_input(step, model)
        out = model.process_step(inp)
        print(f'Step {step}:')
        print({
            '机组排汽流量': out['机组排汽流量'],
            '机组排汽焓值': out['机组排汽焓值'],
            '机组汽轮机效率': out['机组汽轮机效率'],
            '机组冷端总电耗': out['机组冷端总电耗'],
            '机组模型背压仿真误差': out['机组模型背压仿真误差'],
            '机组模型负荷仿真误差': out['机组模型负荷仿真误差'],
            '机组平均风机转速': out['机组平均风机转速'],
        })
        print('-' * 50)


if __name__ == '__main__':
    # main()
    model = Model_Simulation({})
    model.mockData()

